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ABSTRACT 


Two  efficient  algorithms  are  presented  for  obtaining  steady-state 
solutions  of  nonlinear  circuits  and  systems  driven  by  two  or  more 
distinct  frequency  input  signals.  These  algorithms  are  particularly 
useful  in  cases  where  the  steady-state  response  is  either  not  periodic, 
or  is  periodic  but  its  period  is  too  large  for  existing  methods. 

The  first  algorithm  is  applicable  to  any  circuit  or  system  driven 
by  any  number  2  of  input  frequencies.  The  second  algorithm  is 
restricted  only  to  2  input  frequencies  and  is  therefore  significantly 
more  efficient  than  the  first  algorithm.  Both  algorithms  are 
formulated  for  systems  described  by  an  implicit  system  of  nonlinear 
algebraic-differential  equations,  thereby  obviating  the  need  to  write 
state  equations. 

Numerous  examples  have  been  solved  successfully  using  these  two 
algorithms.  A  selection  of  some  of  these  examples  is  given  for 
illustrative  purposes. 
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I.  Introduction 


A  fundamental  problem  in  the  design  of  communication  circuits*  such  as 
modulators  and  mixers,  is  to  calculate  the  steady-state  response  when  the  circuit 
is  driven  by  inputs  having  "P"  distinct  frequency  components  {u^ 
where  P  >  2  [1-2].  For  complete  generality,  we  assume  the  circuit  or  system  is 
described  by  an  implicit  system  of  differential -algebraic  equations  [3]  of  the  form 


fjU.x.yio^t.oigt,. . .u>pt)  =0  ,  j  =  l,2,...,m+n 


(1.1) 


where  x  is  an  n-vector  denoting  the  state  variables,  y  is  an  m-vector  denoting  the 
remaining  non-state  variables  and  f.(*)  contains  p  periodic  input  signals  of  fre- 

J 

quencles  .wg,. • • ,wp,  respectively.  In  Appendix  A,  fj(-)  is  given  by  an  explicit 
formula  which  holds  for  most  circuits  of  practical  interest. 


Is  called  the  steady  state  response  ,  where  the  summation  is  taken  over  a]_l 
possible  frequencies  [5] 


\  ~  mlku>l+m2kul2+- 


..-Hnp^p 


(1.5) 


generated  by  the  frequency  base  to-j ,  *  • .  •  ,wp  • 

Note  that  (1.4)  is  not  an  ordinary  Fourier  series  because  its  frequency  spectrum 
{v-pUg....^}  is  not  harmonically  related.  In  fact,  x$s(t)  is  not  even  periodic  if 
the  frequency  base  {w-j .cjg.- •  •  ,wp}  is  incommensurable  [5].  In  the  mathematical 
literature,  (1.4)  is  called  an  almost  periodic  function. 

Our  objective  in  this  paper  is  to  present  2  efficient  algorithms  for  calculating 
the  steady-state  response  xs$(t) . 
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Current  methods  for  calculating  x$s(t)  can  be  classified  into  4  categories: 

1.  Brute  force  method.  This  approach  solves  (1.1)  by  numerical  integration 
(starting  from  an  arbitrarily  chosen  initial  state  xQ)  until  the  steady  state  is 
reached  [3]. 

Although  this  method  is  quite  general,  it  is  prohibitively  expensive  for 
lightly-damped  circuits  where  it  takes  a  very  long  time  for  the  transient  component 
to  die  out. 

Moreover,  if  the  frequency  base  is  incommensurable,  x$$(t)  is  not  periodic  and 
it  is  difficult  to  determine  when  the  steady  state  has  been  reached. 

2.  Perturbation  method.  This  approach  solves  (1.1)  by  iteration  with  the  initial 
solution  often  chosen  to  be  the  solution  of  a  linearized  equation.  It  includes  the 
Vol terra  series  method  [5-7]  and  the  Picard  iteration  method  [8]. 

Unfortunately,  this  method  works  only  for  almost  linear  circuits  where  the 
nonlinearity  is  often  extremely  weak  (e.g.,  low  distortion  amplifiers).  For  circuits 
which  rely  on  nonlinearity  in  an  essential  way  (e.g.,  modulators  and  mixers)  this 
method  becomes  highly  inaccurate  let  alone  the  fact  that  the  iteration  often  does 
not  converge. 

3.  Harmonic  balance  method.  This  approach  solves  (1.1)  by  approximating  the 
solution  in  a  finite  trigonometric  series  and  then  balancing  all  terms  having  iden¬ 
tical  frequency  components,  often  via  Galerkin's  procedure  [9-10]. 

Although  very  interesting  theoretically,  this  method  is  often  extremely  time- 
consuming  because  the  various  frequency  components  are  estimated  by  multi -dimensional 
Fourier  analysis. 

4.  Shooting  method.  This  approach  solves  (1.1)  by  finding  first  an  initial  state 
Xg(often  via  Newton-Raphon  method)  such  that  the  solution  starting  from  xQ  is 
periodic,  i.e.,  no  transient  component  [11-12]. 

There  are  2  serious  problems  associated  with  this  method. 

(a)  It  can  not  be  used  when  the  solution  is  not  periodic. 

(b)  Even  if  the  solution  is  periodic,  the  period  T  is  often  many  orders  of  magni¬ 
tude  larger  than  the  period  of  the  individual  frequency  components  vk,  thereby 
making  the  numerical  integration  over  this  long  period  T  prohibitively  expensive. 

For  example,  consider 

x$$(t)  =  A.j  cos  o)-|t  +  Ag  cos  ojgt  (1.6) 

The  following  table  lists  several  combinations  of  and  ^  (of  periods  T1  and  T2 
respectively)  which  makes  x  (t)  a  periodic  function.  Also  listed  is  the  period  T 
of  x$s(t)  and  the  ratio  p-j  -  T/T^  and  p2  =  T/T2. 
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Table  1 .  Example  of  and  ^  which  make  x$s(t)  periodic  of  frequency  w. 


cd2(Hz) 

BUH 

v 

c  <x>2 

T-  £<sec) 

1  'l 

A  T 

p2“  T" 

L  l2 

i 

0.23 

6.2832 

27.318 

6.2832(102) 

102 

0.23(102) 

i 

0.233 

6.2832 

26.967 

6.2832(103) 

103 

0.233O03) 

i 

0.2333 

6.2832 

26.932 

6.2832(104) 

104 

0. 2333( 1 04) 

i 

0.2333. -.3 

V  M  J 

6.2832 

26.927 

6.2832(10°) 

10° 

0.233... 3(10°) 

n  digits 

n  digits 

103 

0.233(103) 

6 . 2832 ( 1 0-3) 

0 . 26967 ( 1 0~3) 

6.2832 

103 

0 . 233(1 03) 

104 

0.2333(104) 

6.2832(10"4) 

0.26932(10”4) 

6.2832 

104 

0.2333(104) 

105 

0.23333 ( 1 05 ) 

6.2832(10'5) 

0.26932(10”^) 

6.2832 

105 

0.23333(105) 

10n 

0. 23333.. .3(10n) 

6.2832(10”°) 

0.26927(10”°) 

6.2832 

10° 

0^.233...  3(10") 

1  ^2 

Note  that  when  w-j  s  1,  T  +  «as  n  +  «,  and  when  =  10°,  T-j  -*•  0  but  T  =  6.2832 
as  n  +  «.  Hence,  from  a  numerical  integration  point  of  view,  it  will  take  an 
infinite  amount  of  integration  steps  in  order  to  obtain  the  periodic  solution  xss(t) 
when  n  -*■  «.  Since  p-j  -*-  »  and  p2  -*•  «  as  n  -♦•  <»  in  both  cases,  the  larger  the  values 
of  p.|  and  p2,  the  more  computer  time  will  be  required.  Hence,  p1  and  p2  give  a 
measure  of  numerical  efficiency  of  the  shooting  method.  This  observation  motivates 


the  following: 

Theorem  1 

The  steady  state  response  *ss(t)  in  (1.4)  is  periodic  of  frequency  to  if  each 
frequency  vk  can  be  expressed  as  a  rational  number 
m,. 


vk  =  n 


~  ,  k  =  1,2 . M 


(1.6) 


Moreover,  if  m^  and  nk  are  relatively  prime  integers  for  all  k  «  1,2,...,M, 
then  the  period  T  ^  2ir/w  of  x$$(t)  is  given  explicitly  by 


where* 

n  A  l.C.M.  {n1,n2,...,nM} 
m  ^  G.C.D.  {m^  .n^ . n^} 

*L.C.M.  and  G.C.D.  denote  Least  Common  Multiple  and  Greatest  Common  Divisor, 
respectively. 


(1.7) 

(1.8) 
(1.9) 
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and 


°k  4  "  <i)vk  0  -10) 

where  Tk  =  2ir/vfc. 

Conversely,  if  there  are  at  least  2  frequencies  Vj  and  vk  where  Vj  is  rational 
but  is  irrational ,  then  xss(t)  is  not  periodic. 

Proof.  Rewriting  (1.4)  as 


^jj(t)  —  Xg j  ( V'j  t ,V£t ,  .  .  .  »Vj^t ) 


to  emphasize  the  M  periodic  components  of  frequency  ,v2,. . . ,vM,  we  obtain 
xss(t+T)  =  xss(v1(t+T),v2(t+T),...,vM{t+T)} 

■  *ss  h  (t+2ir(i) >  ^2(1 t+2ir(i)  > » •  •  •  ‘vM(t+2ir(m) 

■  is.  h  <I>  <£>  1  -  €>  > . .  v €» » 

m  t  m  H 


“  5ss(vl(t+  n^OV»Vt+  n^OV  »•  *  *  ’vM*t+  rijjOV) 
3  5ss^vl ^t+NiT] ) »v2^t+N2T2^ ’ * ,vM^t+NMTM^ 


"2,n, 


% 


where 


m, 


V" 


"k  A^S»  -  k  -  1.2 . » 


(1.11) 


(1.12) 


(1.13) 


is  an  integer  in  view  of  (1.8)  and  (1.9).  It  follows  from  (1.13)  and  (1.4)  that 


5ss^t+T)  =  ?ss^lt,v2t,,,*’vMt)  =  Sss^  (1.14) 

Hence  xss(t)  is  periodic  of  period  T.  Moreover,  since  mk  and  nk  are  relatively 
prime,  T  is  the  smallest  period  and  hence,  T  =  2tt/u>. 

Finally,  If  Vj  is  rational  but  vk  is  irrational,  we  can  represent  vk  by  (1.6) 
with  nk  -*•<*>.  Consequently,  n  *  L.C.M.{n^ ,. .  .n^ ,. . .  ,nk». .  .nM>  s  ®  and  x$s(t)  has 
an  infinite  period;  i.e.,  it  is  not  periodic.  ■ 


It  follows  from  Theorem  1  that  if  xs$(t)  is  periodic,  its  frequency  is  given  by 


G. C . D . {m^ ,m2 , . . . ,m^) 

L.C.M. {n^  ,n2>. . .  ,  } 


(1.15) 


and  Its  period  T  is  bounded  by: 

maxtTj .Tg,. .,Tp}  <  T  <  . .nM)  (1.16) 

It  also  follows  from  (1.8)-(1.10)  that  since  m  is  typically  a  small  integer 
(npl  if  all  are  relatively  prime),  the  period  T  can  be  many  orders  of  magnitude 
larger  than  that  of  T^.  Typically,  T  increases  by  an  order  of  magnitude  if  we 
increase  the  number  of  significant  figures  in  representing  the  component  frequencies 
vk,  k  =  1 ,2, . . .  ,p  by  one. 

Observe  that  in  solving  (1.1)  by  numerical  integration,  the  step  size  h  is 
determined  by  the  period  of  the  highest  frequency  component  [3],  namely, 

h  £  g  mi  ntT'j  »T2».*.»Tk}  (1.17) 

It  follows  from  (1.16)  and  (1.17)  that  both  the  brute  force  method  and  the  shooting 

method  are  usually  impractical  when  there  are  multiple  input  frequencies. 

To  overcome  the  problems  associated  with  existing  methods,  we  will  present  two 

new  efficient  algorithms  in  this  paper.  The  basic  idea  in  both  algorithms  is  to 

* 

find  an  initial  state  x(0)  A  xQ  so  that  the  transient  component 

xtr(t)  *  0  for  all  t  >_  0  (1.18) 

regardless  of  whether  the  steady  state  response  x$s(t)  is  periodic  or  not.  In  both 
algorithms,  x*  is  found  by  a  Newton- Raphson  mehtod.  However,  unlike  the  shooting 
method  [11],  (1.3)  is  solved  numerically  only  over  a  small  fraction  of  the  period  T 
(in  the  periodic  case)  per  iteration.  This  is  why  our  algorithms  are  computationally 
quite  efficient. 

The  algorithm  to  be  presented  in  Section  II  is  completely  general  and  is  appli¬ 
cable  regardless  of  the  number  "p"  of  input  frequencies,  provided  p  <  °°. 

The  algorithm  to  be  presented  in  Section  III  is  restricted  only  to  the  2-input 
frequency  case  (p=*2).  We  will  see  that  this  restriction  leads  to  a  significantly 
more  efficient  algorithm  than  that  of  Section  II. 


II.  Almost-Perlodic  Solution  Algorithm  I:  Multiple-Input  Frequencies 

Since  our  algorithm  does  not  depend  on  whether  xss(t)  is  periodic  or  not,  let 
us  assume  that  the  exact  steady-state  response 


(2.1) 
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is  not  periodic  for  the  sake  of  generality.  Consequently,  we  will  call  the  co¬ 
efficients  afc  as  generalized  Fourier  coefficients. 

Even  though  the  number  M  of  frequency  components  in  (2.1)  could  be  quite  large 
(M  may  equal  «!)  in  most  practical  cases. 


A,.  A 


/  ’?2k-l 


0  for  all  k  >  N 


(2.2) 


where  1*1  denotes  Euclidean  norm. 

Hence,  we  will  seek  to  find  an  approximate  solution 


(2.3) 


whose  N  <  M. 

A.  Calculating  ak  when  transient  component  is  zero 

In  section  II-B,  we  will  present  an  algorithm  for  finding  an  initial  state 
Xq  such  that  the  transient  solution  component  *tr(t)  in  (1.3)  is  zero  for  t  >_  0. 
In  this  subsection,  let  us  assume  xQ  has  been  found  so  that  the  solution  of  (1.1) 
starting  from  xQ  is  x(t)  =  xss(t)  for  t  ^  0. 

To  minimize  computation  time,  we  will  often  choose  a  relatively  small  N  so 
that  (2.2)  is  not  necessarily  satisfied.  In  this  case,  the  following  theorem  is 
important: 

Theorem  2.  Properties  of  Generalized  Fourier  Coefficients 
Let  N  <  M  be  an^  positive  integer. 

(a)  For  arrjr  (not  necessarily  optimum)  N-frequency  component  approximation 
N 

xN(t)  £  bQ  +  yb2|M  cos  vkt  +  b2k  sin  vkt) 
to  X  (t)  in  (2.1),  the  mean-sguare  error 

M[xSs(t)-xN(t)]2  A  11m  I  £  lxss(t)-xN(t)l2dt 

is  given  explicitly  by:+ 

+We  define  the  mean  of  x(t)  by 

M{x(t)}  A  lim  i  f  x(t)dt 
T-*»  1  '0 


(2.4) 

I 
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MCXssM-x^t)]2  =  M[xgs(t)]2  -  »aQl2  -  \  l  lakI2  +  la0-b0H2  +  j  l  lak-bkl2  (2.5) 

k-1 

(b)  Among  all  possible  coefficients  {bQ,b1 .  ,b2N}  in  (2.4),  the  coefficients 
which  result  in  a  minimum  mean-square  error  are  precisely  the  first  2N+1  generalized 
Fourier  coefficients;  namely 


bk  =  ak,  k  «  0,1 ,2, . . . ,2N 
(c)  The  minimum  mean-square-error  is  given  by: 

M[x$s(t)-xN(t)]2  -  \  I  latl2 


:n 


O  Z. 

c  k=N+l  ~K 


Proof.  It  suffices  to  prove  the  scalar  case. 

?  N  , 

(a)  M[x$s(t)-xN(t)]^  =  M{x$s(t)  -  bQ  -  £  [b2k_.,cos  vkt+b2|(sin  vkt]r 

k= 1  n 

=  M[xs$(t)]2  -  2b0M[xs$(t)]  -  2M{x$$(t)  £  Cb2k_iCOs  vkt+b2ksin  vkt]} 

N  ^ 

+  b2  +  M{^[b2k_i  cos  ^kt+b2k  sin  vkt]}2 
where  we  have  made  use  of  the  fact  that 


M{sin  vt>  *  0  for  all  v 

M{cos  vt)  3  0  for  all  v  f  0 


(2.6) 


(2.7) 


(2.8) 

(2.9) 

(2.10) 


The  last  term  in  (2.8)  can  be  further  reduced 
N  o  -i  N 

k; 

Substituting  (2.1)  for  x$s(t)  in  the  second  and  third  terms  in  (2.8),  we  obtain 


N  9  ,  N  9  .  2N  7 

^b2k-lC0S  vkt+b2k  sin  vk^}  =  2  ^  (b2k-l+b2k^  3  2  bk  (2J1) 


2b0M[xss(t)]  +  2M{x?s(t)  £  [b2k_1cos  vkt+b2ksin  vkt]} 

N  ksl  2N 

3  2  a0b0  +  k|1(b2k-la2k-1+b2ka2k)  3  2a0b0  +  JjVk 

Substituting  (2.11)  and  (2.12)  into  (2.8),  we  obtain: 

o  9  2N  9  ,  2N  0 

M[xss(t)-xN(t)]  3  M[xss(t)]2  -  2a0b0  -  £  akbk  +  b2  ♦  \  £  b2 

,  9  i  2n  9  ,  ,  2n 

■  M[xss(t)]  -  7  ^  ,  ,J().bo)2  ,  1.  (>k.bk)2 

which  is  just  the  scalar  version  of  (2.5). 
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(2.12) 


(2.13) 


iWIMMk 


(b)  Since  bk  occurs  only  in  the  last  2  terms  of  (2.13),  the  minimum  of  (2.13) 
occurs  when  (2.6)  holds.  M 

(c)  M[x$s(t)]2  =  a2  +  \  l  a2  (2.14) 

Substituting  (2.6)  and  (2.14)  into  (2.13),  we  obtain: 

M[x  (t)-xN(t)]2  =  i  I  a2  (2.15) 

ss  n  ^  k=N+1  K 

which  is  just  the  scalar  version  of  (2.7).  ■ 


Theorem  2  shows  that  regardless  of  the  integer  N,  the  generalized  Fourier  co¬ 
efficients  {aQ,a^ ,. . . .a^nj}  in  (2.3)  can  be  obtained  by  minimizing  the  means-square- 
error  between  xs$(t)  and  xN(t).  Hence,  increasing  the  number  of  frequency  compo¬ 
nents  from  N  to  N+l  does  not  affect  the  previously  calculated  coefficients. 

Since  by  assumption,  xtTlt)  =  0  for  t>0,  we  can  calculate  *ss(t)  by  solving 
(1.2)  numerically.  Assuming  a  uniform  step  size  At,  let  us  calculate  (Z+l)  time 
steps  to  obtain  x  (At),  x.c(2At) ,. . . ,x,c(kAt) ,. . . ,xce ((Z+l )At) ,  where  Z  is  some 
integer  to  be  chosen  later. 

Since  Theorem  2  implies  that  the  ith  components  ak.,  k  =  0,1,. ..,2N,  of  each 
generalized  Fourier  coefficient  ak  can  be  determined  independently  of  the  coef¬ 
ficients  of  the  remaining  components  of  the  vector  xN(t),  it  suffices  for  us  to 
derive  a  formula  for  calculating  these  coefficients  in  the  scalar  case.  To  simplify 
our  notation,  define  the  vectors 


Xss(°) 

xN(o) 

ao 

1 

O 
|<  «TJ 

X$s(At) 

XN( At) 

al 

^(Z) 

xss(z)  a 

xs$(2At) 

>xN(Z)  A 

xN(2At) 

,a  A 

a2 

,  and  a (Z)  a 

a2(Z) 

x$$(ZAt) 

xN(ZAt) 

CM 

fO 

_ _ / 

_a2N(Z)_ 

and  the  Z  *  (2N+1)  matrix 


r  a 


l 

l 

i 


l 

COS  V^At 

cos  2v^ At 


0 

sin  v^At 
sin  2v  At 

•  I 


1  0 
cos  v^At  sin  v^At 

cos  2v^At  sin  2v^At 


1  cos  Zv^At  sin  Zv^At 


cos  Zv^At  sin  Zv^At 


(2.16) 


(2.17) 
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Theorem  3 

The  (2N+1)  generalized  Fourier  coefficients  {aQ,a1 . ,ak>. . . of  xN(t) 
(scalar  version  of  (2.3))  are  given  by 


a  =  a(Z)  +  e(Z) 

where 

a(z)  a  (rTr)_1  rT  xss(z) 


and  e(Z)  is  an  error  vector  satisfying 


(2.18) 

(2.19) 


e(Z)  -*-0  as  Z  (2.20) 

The  (2N+1)  x  (2N+1)  matrix  (F^F)  in  (2.19)  is  non-singular  for  all  positive 
frequencies  »V2**,,,VN  an(*  for  an^  step  size  At  if,  and  only  if 

2nir 


At  / 


lvii'V 


(2.21) 


for  all  i,  k  *  1,2,...,N,  and  for  any  integer  n. 
Proof.  -r 


MOss(t)-xN(t)]2  ■  lim  I 

T-KX) 


Cxss(t)-xN(t)]2  dt 


=  lim 


2-k»  (Z+l  )At 


l  [x  (kAt)-xN(kAt)]2  At, 

l^k*0  SS  N  1 


*  Jo  Cxss(kAt)-xN(kAt)32j  +  e(z) 


(2.22) 


where 

s(Z)  _>  0  denotes  the  error  resulting  from  taking  only  a  finite  number  Z  of  time 
steps. 

Now  substituting  t  =  kAt  in  (2.3)  and  using  the  notations  in  (2.16)  and  (2.17), 
we  can  write 

JQ[xss(kAt)-xN(kAt)]2  =  [x$s(Z)-xN(Z)]TCxss(Z)-xN(Z)] 

=  [xss(Z)-ra(Z)]T[xss(Z)-ra(Z)] 


=  xJs(Z)x$s(Z)  -  aT(Z)rTxss(Z)  -  xJs(Z)ra(Z) 


ss' 

+  aT(Z)rTra(Z) 


-ss' 


(2.23) 


-9- 


(2.24) 


After  adding  and  subtracting  oV^GG^a  «  STG-1a  to  (2.23),  where 
G  A  rTr  and  a  A  rTx$s(Z) 

we  obtain 

kHxss(^t)-xJ,(kAt)]2  *  jJs(2)xss(Z)  -  aT(Z)GG-1o  -  aTGG_1a(Z)  +  aT(Z)Ga(Z) 

+  oV'gG^o  -  oV’a 

=  (a(Z)-G_1o)TG  (a(Z)-G-1o)  +  xJs(Z)x$s(Z)  -  oV’o  (2.25) 

Since  only  the  first  terra  of  (2.25)  depends  on  a(Z),  and  since  G  as  defined 
In  (2.24)  is  clearly  positive  semi -definite,  it  follows  from  (2.22)  and  (2.25) 
that  the  mean-square  error  MCxs.(t)-xf<(t)]^  attains  its  minimum  when 


a(Z) 


G-]a 


(rTr)‘1rTxss(z) 


(2.26) 


It  follows  from  Theorem  2  that  as  Z  ♦  «.  a(Z)  -*•  a  and  hence  e(Z) 

- y  -  - 

The  proof  showing  (r  r) Is  nonsingular  If  and  only  if  (2.21)  holds  involves 

some  cumbersome  determinant  expansions. 


0  in  (2.18). 
jlds  involv 

The  details  are  given  in  Appendix  B. 


Theorem  3  gives  us  an  explicit  formula  for  calculating  (approximately)  the 
(2N+1)  generalized  Fourier  coefficients  {aQ,a1 ,. . . ,a2N)  for  an£  N;  namely. 


a  =  (rTr)'1rT  x$s(z) 


(2.27) 


This  approximate  formula  becomes  exact  as  the  number  Z  of  integration  time 
steps  tends  to  ». 

Since  (2.21)  applies  only  at  a  countable  set  of  isolated  points,  it  is  easy 
to  choose  a  suitable  At  satisfying  (2.21).  Once  At  is  chosen,  r  can  be  calculated 
from  (2.17).  The  vector  £SS(Z)  in  (2.17)  represents  (2N+1)  "samples"  taken  from 
the  exact  steady-state  solution  x  (t)  from  (2.1)  at  a  regular  time  interval  equal 
to  At. 

In  practice,  x  (Z)  is  of  course  calculated  numerically  by  solving  (1.1) 

~5S  *  X 

starting  from  any  initial  state  Xg  which  results  in  a  zero  transient  component. 

Finally,  note  that  (2.27)  gives  the  (2N+1)  generalized  Fourier  coefficients 
of  only  one  component  of  the  state  vector  x$s(t).  Hence,  (2.27)  must  be  applied 

'This  Implicit  system  of  differential -algebra  1c  equationscan  be  solved  using  the 
Backward  Differentiation  Formula  (BOF)  as  described  in  [3]. 


n  times  for  the  "n"  state  variables  in  x.  Since  r  remains  unchanged,  each  calcu¬ 
lation  involves  only  changing  x  (Z)  in  (2.27)  corresponding  to  each  component  of 

~ss(t)* 

B.  Finding  the  initial  state  xQ 

If  we  let  x^  (T^)  denote  the  ith  component  of  x^t)  from  (2.3)  at  any  time 
t  *  Tb ,  then  1 


l.<Tb!  *  a0,  *  |(I1(a2k-lj“svk1b  +  a2k,  s("  ukV 


(2.28) 


when  a^.  denotes  the  ith  component  of  a^.  Substituting  (2.27)  for  a.  in  (2.28), 
we  obtain 

xN.(Tb)«  YT(Tb)ai  «  YT(Tb)C(rTr)“1rTx  (2)D  *  xjs(z)  [r(rTr)_1Y(Tb)], 

i  ■  1,2 . n  (2.29) 

where  x^s  (Z)  denotes  the  xs$(Z)  (as  defined  in  (2.16))  associated  with  the  ith 
component1 of  xs$(t),  and 


y(Tb)  A  [1  cos  slnv^,  ....  cos  vNTb  sin  v^]1 


(2.30) 


Since  r(rTr)_1Y(Tb)  in  (2.29)  remains  unchanged  for  all  i  =  l,2,...,n,  we  can 
combine  all  n  components  of  xN(t)  from  (2.29)  into  a  single  matrix  equation: 


XN^V 

XN0*V 


XN 


XSS1(0)  Xss/At>  xss/2At> 


xec  (ZAt) 


ss 


xss,(0)  xss2(it>  xss2<2m> 


xss  (ZAt) 


xss  (0)  xss  {At)  xss  (2At) 
n  n  n 


x$s  (ZAt) 
n 


?N^Tb*?0  ) 

where  the  (Z+l)-vector 


X(At,Z;  Xq  ) 


’a0(Tb;At,Z,N) 
a-|  (Tb  ;At,Z,N) 
a2(Tb;At,Z,N) 


az(Tb  ;At,Z,N) 


a(Tb;At,Z,N)  (2.31) 


a(Tb;At,z,N)  A  r(rTr)-1 j(Tb) 


(2.32) 


depends  only  on  Tb,At,Z,  and  N  (since  T  in  (2.17)  depends  on  At, Z, and  N)  but  not 


On  Xq  . 
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We  can  recast  (2.31)  into  an  exact  equation  by  introducing  a  slack  variable 
£N(Tb;At,Z);  namely. 


xN(Tb;Xg)  -  X(At,Z;Xg)  a(Tb;At,Z,M)  =  -^(T^-.At.Z) 


(2.33) 


Observe  that  eN(Tb;At,Z)  -*■  0  if  the  following  3  conditions  are  satisfied: 

1 )  Z  =  «. 

2)  x  (t)  is  calculated  by  solving  (1.1)  numerically  with  Xg  as  the  initial 
state,  where  xQ  is  any  initial  state  which  satisfies  (1.18). 

3)  xss(t)  calculated  with  Infinite  precision  and  zero  local  truncation 
error . 

None  of  these  conditions  can  be  exactly  met  in  practice.  Moreover,  even  if 
eN(Tb;At,Z)  =  0,  we  still  can  not  solve  (2.33)  for  xQ  since  we  do  not  as  yet  have 
the  information  (i.e.,  the  coefficients  ag,a^ ,a2,. . .a2N)  needed  in  (2.3)  to  cal¬ 
culate  xN(Tb;xQ  ). 

Instead  of  (2.33),  however,  suppose  we  define  the  following  system  of  "n" 


where  x.  (t*,Xg)  denotes  the  1th  component  of  the  complete  solution  x(t)  =  xtr(t) 

+  x$$(t)  0*1)  starting  from  the  initial  state  Xg,  for  t  =  0,  At,2At,. . . ,ZAt, 

and  T^.  Observe  that  unlike  (2.31),  both  x(Tb;Xg)  and  X(At,Z;xQ)  can  be  calculated 
by  solving  (1.1)  numerically. 

Now  at  Xg  =  x0  ,  xtr(t)  =  0  for  t  _>  0  (by  definition)  and  hence  we  can  write 


*<Tb  ;?0  } 


?ss^Tb’?0  ^  *  ?N^Tb;?0  ^  +  ^M-N^b  ;?0  ^ 


(2.35) 
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where  ?ss^b;~0*^  =  Sss^V  and  ^N<Tb^O*)  =  ?N^Tb^  as  defined  in  (2-1)  and  (2-3) 
respectively;  and  where  xM_N(Tb ;xQ  )  denotes  the  remaining  (M-N)  terms  of  xss(Tb) 
which  have  been  excluded  from  x^T^-  Consequently,  at  =  Xq  ,  (2.34)  can  be 
written  as  follows: 

Xjj(l"b»^o  ^  ~  X(At,Z;Xg  =  ”?M-N^b’?0  ^  (2.36) 

Comparing  (2.33)  and  (2.36),  we  find 

~N^b’^’^  =  ?M-N^b’?0  ^  (2.37) 

Equation  (2.37)  is  remarkable  because  it  says  that  eN(Tb;At,Z)  0  when 
N  -*•  M.  In  other  words,  if  the  exact  steady  state  response  x$s(t)  in  (2.1)  has 
only  M  <  »  frequency  components,  and  if  we  choose  N  =  M  in  (2.3),  then 
£N(Tb;At,Z)  »  0  and  (2.31)  becomes  exact  for  any  Z. 

Indeed,  when  M  =  N  and  Z  *  2M  +  1 ,  r  becomes  a  square  matrix  and  the  gen¬ 
eralized  Fourier  coefficients  can  be  calculated  exactly  from  (2.27): 

a  =  r_1xss(Z)  (2.38) 

Similarly,  (2.32)  in  this  case  (Z=2M+1)  reduces  to 

a(Tb;At,Z)  -  rTY(Tfa)  (2.39) 

Of  course  in  practice,  we  will  normally  choose  N  «  M  in  order  to  save  compu¬ 
tation  time.  This  choice  is  often  necessary  anyway  because  M  *  °°  for  most  prac¬ 
tical  circuits.  Fortunately,  the  amplitudes  of  the  higher-order  terms  [5]  in 
many  practical  circuits  satisfy  (2.2)  so  that  the  error  vector  e(Tb;At,Z)  remains 
relatively  small  even  though  N  «  M. 

Let  us  summarize  the  preceding  observations  as  follow: 

Remarks: 

A 

1.  The  solution  Xq  of  the  nonlinear  equation 

F(x0;Tb,At,Z,N)  *  0  (2.40) 

as  defined  in  (2.34)  for  fixed  At  and  Z  represents,  a  good  approximation  to  x^ 
provided  the  number  of  frequency  components  N  and/or  the  number  of  time- step 
samples  Z  are  chosen  to  be  sufficiently  large.  In  particular, 
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x0  -*■  xg  as  N  -*■  “  and/or  !■*■<*> 


(2.41) 


2.  Since  (2.40)  is  not  given  in  closed  analytical  form,  it  must  be  caluclated 
numerically  for  each  Xg,  At,  and  Z. 

3.  Equation  (2.40)  can  be  solved  for  xQ  (for  fixed  At  and  Z)  by  the  Newton-Raphson 
method  [3]:  . 


?0J+1)  =  ?0j)  ‘  ^F(?0j))]’1  f(?Oj);Tb’At,Z’N) 
iil  ^^(xn’^h *^t»Z ,N) 

^iA3)n^4- -  ,d) 

~u  CO  CO 


(2.42) 


(2.43) 


denotes  the  Jacobian  matrix  of  F(xg;Tb,At,Z,N)  at  xQ  »  XgJ^.  This  can  be  evalu¬ 
ated  by  the  method  given  in  Section  II-C. 

★  ★ 

4.  Once  the  initial  state  Xg  is  found,  we  solve  (1.1)  numerically  with  Xg  as 

the  initial  state  to  obtain  x  (At),  x  (2At),...,  x  (ZAt).  Substituting  the  ith 

•*jj  —55  —55 

component,  i  *  1,2 . n,  of  these  data  into  (2.27),  we  obtain  the  first  (2N+1) 

generalized  Fourier  coeeflcients  aQ  ,a1  ,...,a2N  of  the  ith  component  x$s  (t)  of 

the  steady  state  response  x$s(t).+  ’lWsteady-state  response  xN(t)  at  anytime 

t  *  T,  can  now  be  obtained  by  calculating  (2.3)  at  t  »  T.. 

J  (  A  \  J 

C.  Evaluating  the  Jacobian  matrix  Jr(x^Jj) 

Since  the  most  tlme-consuminq  part  In  solving  for  xn  via  the  Newton-Raphson 

-U  /  2\ 

method  is  the  numerical  calculation  of  the  Jacobian  matrix  JpUg  )»  it  is 
essential  to  develop  efficient  computational  methods.  Taking  the  Jacobian  of 
F(x0;At,Z,N)  in  (2.34),  we  obtain 


JF{?0 } 


9;(Tb ;x0) 


X  *x(j) 
•CO  ~0 


Z  3x(kAt;xn) 

L  °k<v^>  - 


x  -x<J> 
'CO  -0 


(2.44) 


Hence,  we  need  to  calculate 


3x(t;Xg) 

3xn 


i?o3*6 


at  t  *  0,At,2At,...,ZAt,  and  Tb 


(2.45) 


These  (Z+2)  n  x  n  matrices  can  be  calculated  by  the  numerical  differentiation 
method  described  in  Section  17-5-2  of  [3].  If  the  circuit  associated  with  (1.1) 

If  xss(t)  Is  periodic  and  Its  period  T  is  not  too  large,  we  can  replace  this  step 
by  numerically  solving  (1.1)  from  t  *  0  to  t  *  T  with  xn*  as  the  initial  state. 
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is  given,  the  most  efficient  method  for  calculating  (2.45)  is  the  sensitivity 
network  approach  given  in  Section  17-5-3  of  [3], 

However,  if  (1.1)  is  available  only  analytically*,  the  sensitivity  network 
approach  is  not  applicable.  In  this  case,  the  following  method  is  much  more 
efficient  and  accurate  than  that  of  numerical  differentiation: 

In  vector  form,  (1.1)  becomes 


f  (x ,x ,y t ,(^>2^  t  *  ■  *  *u)pt )  =  0 


Applying 
where  fac 

x0  =  x0  ’  we  obtain 


g  Taylor  Expansion  about  (x^(t).y^(t))  at  the  jth  stage  of 
X^(t),y^(t))  denotes  the  solution  of  (2.41)  with  initial  s 

^  iiia  nht1,  i  n 


(2.46) 

the  iteration, 
state 


f(x^^(t),x^^(t),y^^(t),£o1t,oj2t,...,ajpt)  + 


+  0(l^(t)l2,Mt)l2,HY(t)S2)  =  0 


3f  3f  3f 
3X  3x  3y 


n(t) 

n(t) 

y(t) 


(2.47) 


where 

n(t)  A  x(t)  -  x^(t),  y(t)  A  y(t)  -  y^(t)  (2.48) 

The  first  term  in  (2.47)  is  Identically  zero  because  (x^(t),y^(t))  is  a  solu¬ 
tion  of  (2.46).  Neglecting  the  higher-order  terms,  (2.47)  can  be  recast  as 
follow: 


r ,  — 1 

-1 

n(t) 

3f  3f 

r»!i 

1 

>-l 

_ 1 

3x  3y 

i 

•  Sfl- 

i _ 

n(t)  a 

A(j)(t) 

B(j)(t) 

(x(j)(t),y(j)(t)) 

- 

n(t) 


(2.49) 


The  first  component  equation  of  (2.49)  is  a  linear  time-varying  differential 
equation 


+0ur  algorithms  in  this  paper  are  valid  for  an£  equation  of  the  form  (1.1),  which 
need  not  be  associated  with  a  circuit. 
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r 


n(t)  *  A^(t)  n(t) 

where  A^(t)  is  an  n  x  ^-matrix \funotton~of  time.  We^ilX  henceforth  refer  to 
(2.50)  as  the  variational  equation  associated  with  (2.46). 

The  solution  of  (2.50)  corresponding  to  any  initial  state  n(Q)  is  given  by 


(2.50) 


[13] 


n(t)  =  $^(t)  n(0) 


(2.51) 


where  4>^(t)  is  the  fundamental  matrix  solution  of  (2.50).+  If  we  choose 


n(0)  *  [0  0  ...  0  nk(0)  00  ...  0]T 

then 

n^t)  *  *jjJ(t)  nk(0),  j  =  l,2,...,n 


(2.52) 

(2.53) 


where  $..(t)  denotes  the  ikth  element  of  $(t).  It  follows  from  (2.52)  that 

1 K  ~ 

3x.(t)  n^t) 

=  "  *ik(t) 

Hence  we  have  proved  that 


(2.54) 


3x(t;xQ) 

3?0 


x  =x(j) 
00  00 


=  4>(j)(t) 


(2.55) 


It  follows  from  (2.54)  that  0p(x^)  in  (2.44)  can  be  calculated  accurately 
in  3  steps: 

1)  Form  the  variational  equation  (2.50)  at  each  iteration. 

2)  Calculate  the  fundamental  matrix  solution  $^(t)  of  (2.50). 

3)  Calculate  (2.44). 


D.  Initialization  Guidelines 

To  initiate  the  algorithm  for  finding  the  Initial  state  xQ  ,  it  is  necessary 
to  choose  the  5  parameters  N,  Z,  At,  Tb  and  xQ  for  constructing  the  nonlinear 
equation  (2.34).  Since  a  good  choice  of  these  parameters  depends  on  both  the 
nature  of  the  problem  [number  of  state  variables,  degree  of  nonlinearity, 
amplitudes  of  input  signals,  number  of  Input  frequencies,  etc.)  and  the  computer 
being  used  (word  length,  single  or  double  precision,  etc.),  we  can  only  offer 

+The  jth  column  of  the  fundamental  matrix  solution  Is  simply  the  solution  of  (2.51) 
with  the  initial  state 

n  (0)  *  [0  0  ...  0  1  0  ...  0]T 

^ —  jth  position 
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some  guidelines  which  have  been  found  useful  in  our  numerous  numerical  experiments 
conducted  using  our  algorithm. 

a)  Choice  of  N. 

Recall  N  £  M  Is  the  number  of  frequency  components  used  in  the  truncated 
steady  state  solution  xN(t)  In  (2.3).  For  typical  communication  circuits 
(amplifiers,  mixers,  modulators,  etc.)  the  number  of  significant  frequency  compo¬ 
nents  is  usually  known  from  previous  analysis  and  N  should  be  chosen  to  include 
all  such  components. 

If  the  number  of  significant  frequency  components  is  not  known  from  previous 
experience,  we  simply  make  an  intuitive  guess.  If  this  guess  is  unrealistic,  it 
will  show  up  in  the  subsequent  error  analysis  (to  be  discussed  in  Section  II-E) 
and  we  will  have  to  repeat  the  analysis  with  a  larger  N. 

b)  Choice  of  Z,  At,  and  Tk. 

Recall  that  At  is  the  uniform  sampling  step  size  used  in  "sampling"  the 
numerical  solution  of  (1.1)  and  Z  is  the  total  number  of  samples  taken.  It  is 
important  to  note  that  "At"  is  not  the  same  as  the  integration  step  size  "h" 
used  in  solving  (1.1). 

In  most  of  our  numerical  experiments,  we  solve  (1.1)  using  a  4^-6^)  order 
BDF  algorithm  [3]  with  a  step  size  h  =  Tmin/50,  where  Tm-n  is  the  smallest  period 
of  the  N  frequency  components.  This  choice  usually  gives  a  very  accurate  numerical 
solution  for  x(t). 

Our  sampling  step  size  At  is  usually  chosen  within  the  range 


7h  <  At  <  13h 


(2.56) 


provided  (2.21)  is  satisfied.  In  practice,  ill-conditioning  could  occur  if  At  is 
chosen  to  be  too  small,  or  if  it  contains  some  frequency-  components  vj  and  such 
that  f Vj  -  ( =  0.  (See  Appendix  B) 

Although  Theorem  3  shows  that  the  generalized  Fourier  coefficients  can  be 
calculated  exactly  only  If  Z  -*■  «  (see  Eq.  (2.20)),  our  numerical  experiments  show 
that  good  results  can  be  obtained  in  many  practical  cases  with  a  considerably 
smaller  Z.  In  particular  we  have  found  the  following  range  to  be  adequate  for 
the  many  examples  we  have  tried  so  far: 


(2N+1 )  <  Z  <  2(2N+1) 


(2.57) 
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Finally,  the  choice  of  is  somewhat  arbitrary  as  it  does  not  affect  the 
theory  in  Sec.  II-A  from  which  our  algorithm  is  based.  However,  since  the  gener¬ 
alized  Fourier  coefficients  are  estimated  by  samples  taken  over  the  time  interval 
[O.ZAt],  we  choose 


Tb  >  ZAt 


(2.58) 


so  that  the  data  x(Tb;Xg)  would  not  be  redundant, 
c)  Choice  of  Xg. 

To  assure  and  to  hasten  the  convergence  of  the  Newton-Raphson  iteration,  it 
is  desirable  to  pick  a  good  Initial  guess  Xg^.  Unfortunately,  no  intuitive  guide¬ 
lines  are  available  especially  when  the  steady  state  solution  xss(t)  is  not 
periodic. 

One  approach  which  has  worked  well  for  our  examples  is  to  replace  the  input 

A  A  A 

frequencies  ,v2,. . . ,\>N>  by  an  approximate  set  of  frequencies  (v1 ,v2,. . . ,vN)  so 
that  the  associated  steady-state  waveform  is  periodic  with  a  relatively  small 
period  T  =  27r(n/m),  where  m  and  n  are  defined  in  (1.8)-(1.9)  and  is  bounded  by 
(1.16).  Using  this  approximate  set  of  frequencies,  we  then  apply  the  shooting 

method  [11],  or  any  otner  efficient  method  for  findinq  x„  for  periodic  solutions, 

*  *  f  o') 

to  calculate  xQ  .  We  then  take  this  approximate  xQ  as  our  initial  guess  Xg  . 

If  we  let  in  „  A  max{m1,m,,...,mM}  and  n.  A  min{n, ,n0 . nu},  then  (1.7) 

max  —  i  c  N  min  —  id  N 

suggests  the  following  algorithm  for  reducing  T: 

(1)  If  v,.  *  m,.  is  an  integer  for  all  k  =  1.2.....N,  then  we  increase  m _ 

K  K  a  iuoXa 

until  it  is  not  a  prime  number  and  then  increase  v.,  k  =  1,2,...,N,  until  m  /v. 

I\  HKS  A 

is  an  integer. 

Example  1.  let  ^  *  2,  =  3,  and  v3  =  7.  Then  {n^.mg.nig}  =  {2,3,7}  and  we 

have  from  (1.7) 


T  =  mTT2T3TTT  *  2ir 

A  A 

Now  increase  n#  „  =  7  to  v,  *  m,  =  8,  and  then  increase  v,  to  4. 
max  A  A  o  o  c 

associated  with  (v-j  .Vg.v-j)  *  (2,4,8)  is 

T  *  5.C.D’.T2,4,8)  = 


The  new  period 


(2.59) 


(2)  If  vk  *  rn^/n^  is  not  an  integer,  we  first  change  mk  as  in  (1)  and  then 
change  nk,  k  *  1,2,...,N  until  it  becomes  a  multiple  of  nmin. 

Example  2.  Let  *  2/5,  Vg  *  3/8,  and  v-j  a  7/9.  Then  {m^m^.mj}  -  {2,3,7} 
as  in  Example  1  and  {n^.ng.n^}  *  {5,8,9}.  From  (1.7),  we  find 
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T  =  2tt 


L.C.M.{5.8.9} 

G.C.D.{2,3,7} 


=  720ir 


/V  A 

Since  nmin  =  5,  we  change  n 2  and  n3  to  n2  -  10,  and  n3  =  10  so  that 


(2.60) 


T  = 


2ir 


L.C.M. {5,10,10} 
G.C.D. {2,4,8} 


IOtt 


(2.61) 


Note  that  dramatic  reduction  in  period  from  720ir  to  10it!  1 

(3)  If  is  an  irrational  number,  we  first  approximate  it  by  a  rational 
number  and  then  proceed  as  in  (2). 

Example  3.  Let  =  0.404040...,  v2  -  0.375010101...,  and  v3  =  0.7777...  . 

We  can  approximate  v.j,v2,  and  v3  fay  =  2/5,  v2  =  3/8,  and  \> 3  =  7/9  and  then  pro¬ 

ceed  as  in  Example  2.  Note  the  period  changes  from  T  =  »  to  T  *  lOir. 

E.  Termi nati on  Gui del i nes 

Since  our  choice  of  N  may  not  be  realistic  in  the  sense  that  one  or  more 

significant  frequency  components  may  have  been  inadvertently  excluded  from  (2.3), 

our  algorithm  does  not  terminate  when  the  Newton-Raphson  iteration  in  (2.42) 

* 

converges  to  an  initial  state  xQ  .  We  must  further  validate  our  answer  as 
follows: 

(1)  If  the  steady-state  response  x„(t)  is  periodic  with  a  reasonably 

*  * 

small  period  T,  we  simply  solve  (1.1)  numerically  for  x(t,xQ  )  (with  xQ  as 
initial  state)  from  t  =  0  to  t  *  T  and  verify  that  x(0;xQ)  =  x(T;Xq) 

C2)  If  the  steady-state  response  xs$(t)  is  not  periodic,  or  if  it  is 
periodic  with  an  unreasonably  large  period  T,  we  can  carry  out  the  following 
heuristic  validation  procedure  in  view  of  (2.7)  of  Theorem  2. 

(a)  Solve  (1.1)  numerically  for  x(t,XQ  )  (with  Xq  as  initial  state)  from 
t  =  0  to  t  *  T,  where  T  is  chosen  to  be  sufficiently  large. 

(b)  Solve  for  xN(kAt,xQ)  using  (2.31)where  Tb=  kAtand  theentries  x$s  (kAt)  in  the 
matrix  X(At,Z;xQ*)  are  substituted  by  x^. (t;xQ  ),  t  =  0,  At, 

(c)  Calculate  the  error 


1 4{  Lxj  (kAt; ? o* ) “ *N .  (k At  i xQ* ) ] 


(2.62) 


4»  * 

Jie  could  reduce  T  further  by  decreasing,,  (instead  of  increasing)  n2  and  n3  to 

n2  *  5  and  n3  =  5.  However,  the  v2  and  v3  no  longer  represent  a  good  approximation. 
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for  each  component  j  =  1 ,2,. . . ,n. 

If  maxfe-j .eg** •  •  »£n)  is  smaller  than  some  perscribed  tolerance,  stop.  Other¬ 
wise,  increase  N  and/or  Z  and  start  all  over  again. 

F.  Summary  of  Multi -Frequency  Algorithm 

Step  0.  Specify  the  4  parameters  N,  Z,  At,  and  Tb  (See  Section  II-D)  and 
calculate  the  vector  a(Tb;At,Z,N)  using  (2.32). 

Set  j  =  0. 

Step  1.  Choose  initial  state  Xg  =  x^  (for  j  =  0,  see  Section  II-D)  and 
solve  (1.1)  numerically  to  obtain  x(At),x(2At),...,x(ZAt). 

Calculate  F(xg-,Tb,At,Z,N)  from  (2.34).  If  llF(xg;T,At,Z,N)l  <  e 
whose  e  is  a  sufficiently  small  positive  number,  call  xQ  =  xQ  and 
go  to  Step  4. 

Step  2.  Compute  the  Jacobian  matrix  Jp(x^)  in  (2.43).  (See  Section  II-C). 

Step  3.  Compute  Xg'**1^  via  the  Newton-Raphson  iteration  (2.42). 

Go  to  Step  1  with  j  -►  j  +  1 . 

•ft 

Step  4.  Solve  (1.1)  for  x(t)  with  xQ  as  initial  state  from  t  =  0  to  t  =  Tb 
where  Th*  period  if  xss(t)  is  periodic,  orTb  is  a  sufficiently  large 
number. 

Case  1.  x$s(t)  is  periodic  with  small  period  T: 

Calculate  =  x^O.Xg*)  -  Xj(T;xQ  ),  j  =  1,2 . n 

Case  2.  xss(t)  is  not  periodic  or  is  periodic  with  large  T: 

Calculate  e^.  using  (2.62),  j  =  l,2,...,n. 

If  max{e1 ,e2,. . . ,en}  >  eQ 

where  eQ  is  a  sufficiently  small  positive  number,  increase  N  and/or 
Z  and  repeat  Steps  0-4. 

Step  5.  STOP. 


G.  Illustrative  Examples 

Numerous  examples  have  been  solved  successfully  using  the  2  algorithms  pre¬ 
sented  in  Sections  II  and  III.  Because  of  its  widespread  interest,  let  us  apply 
the  preceding  algorithm  to  solve  the  forced  Duffing's  equation  [14-15]: 


(2 


This  equation  arises  in  many  physical  problems  (e.g.,  ferro- resonance  circuits) 
and  is  known  to  exhibit  many  interesting  phenomena;  including  subharmonic,  almost 


.63) 
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periodic,  and  chaotic  solutions  [16]. 

To  apply  our  algorithm,  let  us  recast  (2.63)  into  the  form  of  (1.1),  which 
in  this  case  is  just  the  state  equation 


X1  =  x2 

*  3 

x2  =  -kx2  -  c^  -  c2x^  +  f(t) 


To  be  specific,  let  us  choose  a  3-frequency-input  signal 


(2.64) 


f(t)  =  A-j  cos  u^t  +  A2  cos  w2t  +  A3  cos  U)3t  (2.65) 

and  k  =  0.1,  c-j  =  2.0  and  c2  =  1.0. 

We  have  solved  (2.64)  using  many  different  combinations  of  amplitudes  and 
frequencies,  4  of  which  are  listed  in  Table  2. 

Table  2.  Four  combinations  of  A.,  and  w.  and  their  respective  periods. 


case 

Ai 

A2 

A3 

wi 

w2 

“3 

B1 

m 

B9 

■ 

0.4 

l 

0.35 

0.155 

6.283 

17.951 

40.54 

400tt 

H 

0.4 

0.4 

0.4 

l 

0.85 

0.170 

6.283 

7.392 

36.96 

200tr 

H 

0.5 

0.5 

KH 

l 

0.35 

0.155 

6.283 

17.951 

40.53 

400tr 

■ 

■ 

■ 

B 

l 

0.85 

0.170 

6.283 

7.392 

36.96 

200tt 

From  previous  experience,  we  know  all  frequency  components 
vk  =  mlka)l  +  m2kw2  +  m3ka>3 

with 

lmlkl  +  lm2k '  +  fm3k*  -  3 

are  likely  to  be  non-negligible.  Since  these  are  30  frequency  components  satis 

fying  (2.66),  we  choose  N  =  30  in  (2.3).  Applying  the  preceding  algorithm  with 

★ 

At  *  11  T.j/50,  Z  *  1.5(60)  =  90,  and  Tb  =  23  T^ ,  we  obtain  the  initial  state  xQ 
listed  in  Table  3  corresponding  to  the  4  cases  in  Table  2.  Also  listed  is  the 
error  ej  calculated  using  (2.62) 


(2.65) 

(2.66) 
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Table  3.  Initial  state  computed  using  multi -frequency  algorithm  with  a  6th  order 
BDF  algorithm  [3]. 


xQ*  A  x*(0) 

Error 

case 

x/(0) 

x2*(0) 

el 

e2 

1 

0.69667 

-0.18304 

0 . 53(1 0-3) 

0.99(10'3) 

2 

0.78298 

-0.13834 

0.14(1 0”^ ) 

0 . 22 ( 1 0~2 ) 

3 

0.82931 

-0.32269 

0.92(1Q“3) 

O 

O 

l\> 

4 

0.81562 

-0.46932 

0.33(10'2) 

0.57(10-2) 

Using  the  initial  states  from  Table  3  and  (2.27),  we  have  calculated  the  60 
generalized  Fourier  coefficients  a1  ,a2... .  ,ag0  in  (2.3)  corresponding  to  N  =  30  for 
cases  1  and  2.  The  waveforms  of  x^t)  for  these  2  cases  are  plotted  (using  (2.3))  as 
the  solid  waveforms  in  Figs.  1(a)  and  2(a),  respectively.  As  a  check  over  the 
accuracy  of  our  solutions,  we  solve  (1.1)  using  the  same  initial  states  and  the 
solution  at  each  integration  time  step  is  shown  as  "dots"  in  Figs.  1(a)  and  2(a) 
respectively.  Mote  the  remarkable  accuracy  in  both  cases. 

To  compare  the  amplitudes  of  the  30  frequency  components,  we  use  (2.27)  to 
plot  the  discrete  frequency  spectrum  for  these  2  cases  in  Figs.  1(b)  and  2(b), 
respectively. 

Finally,  to  obtain  a  measure  of  the  rate  of  convergence  of  the  Newton-Raphson 
iteration  (2.42),  the  error 


■«»*/ 


Fl(?0 


{j);Tb  ,At,Z,N)+F^(x£j);Tb  ,At,Z,N) 


(2.67) 


at  each  iteration  is  plotted  in  Fig.  3  for  cases  1  and  2,  respecti vely.  Mote  that 
both  converges  rapidly  in  4  iterations. 


III.  Almost-Periodic  Solution  Algorithm  2:  Two  Input  Frequencies 

In  this  section,  we  assume  the  circuit  or  system  is  driven  by  no  more  than  2 
frequencies;  i.e.,  P  <_  2  In  (1.1).  Hence,  let  us  rewrite  (1.1)  and  (1.5)  as 
f ol 1 ows : 


f(x,x,y;Wlt,u2t)  a  0  (3.1) 

~  ^  rn2ktJL>2 9  ^  s  ^  *^  *  *  •  •  (3.2) 
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Substituting  (3.2)  into  (1.4)  and  making  use  of  standard  trigonometric  identities, 
we  can  recast  the  steady-state  response  x$s(t)  a$  follow: 


-ss 


=  ?0  +  j^?2k-l  cos^mlk“l4m2k“2^t  +  ?2k  sin^mlku*l+m2ku,2)tJ 

M  r 

=  aQ  +  J  a2k.-,  [(cos  mlkw1t)(cos  m2k(jj2t)-(sin  m-j^tKsin  n^oigt)] 
+  a2k[(sin  m^kai1t)(cos  m2ko)2t)+(cos  m^k^t)(sin  m2kuj2t)]^ 

M  r 

=  '°  +  k-1  I  ^~2k_1  C°S  mlkU3lt+~2k  sin  mlku>lt-*  cos  m2kw2t 


+  [a2k  cos  m]|<(JJ-jt"a2k-l  sin  mikult^  sin  m2k‘°2tJ 


If  we  let  B  denote  an  integer  bound  such  that 

I ml  k I  +  l^kl  -  B 


(3.3) 


(3.4) 


then  the  number  M  of  non-zero  frequency  components  vk  is  given  in  Table  4  for 
B  38  1 ,2, ...  ,10. 

Table  4.  The  integers  M,  2M  +  1 ,  and  2B  +  1  as  a  function  of  B. 


B 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

M 

2 

6 

12 

20 

30 

42 

56 

- 

72 

90 

no 

2M+1 

5 

13 

25 

41 

61 

85 

113 

145 

181 

221 

2B+1 

3 

5 

7 

9 

11 

13 

15 

_ 

17 

19 

21 

This  table  can  be  easily  verified  hy  counting  the  number  of  solid  dots  subtended 
by  an  isosceles  triangle  of  base  length  2B  in  Fig.  4.  For  example,  we  can  enum¬ 
erate  the  following  frequency  components  when  B  =  3: 

v-j  —  u)-j ,  v2  —  u^2 ,  —  2oj -j ,  —  2oj2,  Vg  —  3cj-| ,  Vg  —  3 to2»  Vy  —  ^  ^2» 

Vg  —  (ii^  —  Ug,  Vg  —  (ii^  ^  2u2,  v^q  *  -  2(^2 1  V|  ^  —  2a)-j  ^  Wg  >  and  ~  2w^  —  u)2* 

Hence,  M(3)  *  12. 

Observe  that  all  solid  dots  on  the  m^-axis  in  Fig.  4  denote  harmonics  of 
w.|.  Likewise,  those  on  the  m2k-axis  denote  harmonics  of  to2.  All  other  solid  dots 
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denote  intermodulation  components.  In  particular,  all  solid  dots  on  a  horizontal 
line  through  m2k=N,  N=  1,2,...,  correspond  to  frequency  components  of  the  form 
vk  -  +  Na^-  Hence,  if  we  regroup  all  frequency  components  in  (3.3) 

corresponding  to  dots  on  a  horizontal  line  together,  we  can  recast  (3.3)  into 
the  form 

(3.5) 

where  gQ(t),  g-j(t),  g2(t) ,. . . ,g2B(t)  contain  only  cosine  and  since  components 
which  are  harmonics  of  oj-j  and  are  therefore  all  periodic  functions  of  period 
Tj  =  Since  this  observation  is  the  basis  of  Algorithm  2,  we  will  restate 

it  as  a  theorem: 

Theorem  4. 

The  steady-state  response  x$s(t)  in  (3.3)  which  contains  (2M+1)  generalized 
Fourier  coefficients  can  be  recast  into  the  form  of  (3.5)  containing  only  (2B+1) 
coefficient  functions  of  time  gQ(t),  g-j(t),  g2(t) ,. . . ,g2B(t)  which  are  all 
periodic  of  period  T1  =  2-n/u^ . 

A  comparison  between  the  number  of  coefficients  describing  (3.3)  and  (3.5)  is 
given  in  Table  4.  Observe  that  (3.5)  has  much  fewer  coefficients  compared  to  that 
of  (3.3)  specially  for  large  B.  For  example,  when  B  =  8,  Algorithm  1  from  Section 
_II  would  entail  solving  for  145  generalized  Fourier  coefficients,  whereas  only 
17  coefficient  functions  need  be  specified  in  (3.5).  Our  objective  in  this 
section  is  to  develop  a  new  algorithm  which  takes  full  advantage  of  this  remark¬ 
ably  concise  form  of  solution. 

A.  Calculating  qk(0)  when  transient  component  is  zero. 

In  Section  III-B,  we  will  present  an  algorithm  for  finding  an  initial  state 
Xg  such  that  the  transient  component  *tr(t)  in  (1*3)  is  zero  for  t  >  0.  In  this 
subsection,  let  us  assume  Xg  has  been  found  so  that  the  solution  of  (3.1) 
starting  from  xQ  is  x(t)  *  x$s(t)  for  t  >  0,  where  x$s(t)  is  given  by  (3.5). 

For  reasons  that  will  be  clear  in  Section  III-B,  we  need  to  derive  a  rela¬ 
tionship  for  calculating  gk(0),  k  =  0,1 ,2,. . . ,2B,  in  terms  of  "(2B+1)"  samples 
xss(0),  ?ss^V*  *ss^2T1  )»••••  ^s^V  taken  at  Tj  =  2-rr/oj-,  intervals.  Since  each 
component  x„  (t)  of  xec(t),  1  =  l,2,...,n,  can  be  calculated  separately,  it 
suffices  for  us  to  derive  the  1th  component  g^  k(0)  of  gk(0). 

Substituting  t  ?  0,  T-j , 2T-j , . . .  ,2BT^  into  (3.5) 
and  using 

9i  ,k^Ti )  a  9-|^(0)*  k  =  1 .2, . . .  ,2B  (3,6) 

in  view  of  Theorem  4,  we  obtain 
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xss,(0!  ■  91,0<°>  *  J,  91.2k-l<°> 

xSSf(V  ’  91,0(0)  +  j,  f9i,2k-l(0)  cos  k“2Tl  +  9i,2k(0)  s'"  k“2Tl) 

Xss/2T,)  ’  91,0<°>  *  j,{9f,2k-l(0)  cos  2k“2Tl  *  9i,2k(0)  sin  2k“2Tl) 
•  •  • 

xss,‘2BTl>  =  91,0(0)  ♦>  j,{s,t2k.,«»  “*»8“2V  9j>2k(0)  sm2kB«2T^ 

Equation  (3.7)  consists  of  2B+1  equations  in  terms  of  the  2B+1  coefficients 

9i ,0<°> *9i  ,1 ^  ,9i ,2^  »• ‘ * *9i  ,2B^  * 

If  we  define  the  (2B+1 )- vectors 


xss.(°) 

9i 

XS5,<V 

9i 

,1<°> 

;ss,<B)  a 

xSs/2Tl> 

and  gn  A 
~Bi 

9i 

O 

C\J  •  • 
* 

xss,<2BV 

9i 

,2B(0) 

and  the  (2B+1)  x  (2B+1)  square  matrix 


1  1 


ft(B)  A 


1  cos(^2T])  sin(o^T-j) 

1  cos(2oj2T1  )  sin(2u^T1) 


1  0 
cos(Bu^T-| )  sinfB^Tj) 

cosUBu^T.])  sin(2Bo!2T1) 


2  2 

1  cosfBu^Tj)  slnfBugTj)  •••  cos(B  u^)  sin(B  u^T-)) 


1  cosUBwgl^)  sin(2Ba)2Ti) 
then  (3.7)  assumes  the  condensed  form 


cosUB^T-j )  sin(2B2o)2T-|) 


fl(B)gR  *  x  (B),  i  *  1,2,. ..,n 
~  ~Bi  ~ssi 


(3.7) 


(3.8) 


(3.9) 


(3.10) 


-25- 


Theorem  5 


The  2B+1  coefficients  g^  Q(0) ,g^  1(0),gi  2(0) , . . 2B(0)  describing  the 
steady-state  response  (3.5)  can  be  calculated  exactly  from 


i  =  1 ,2,. . . ,n 


(3.11) 


The  matrix  Q(B)  is  non-singular  if,  and  only  if,  there  does  not  exist  an  integer 


L 


2 


such  that 

p-  .  L,  «  1,2,. ...2B 
"l  L1  1 


(3.12) 


Proof.  Eq.  (3.11)  follows  directly  from  (3.10).  The  proof  that  (3.12)  is  a 
necessary  and  sufficient  condition  for  n(B)  to  be  non-singular  is  given  in 
Appendi x  C .  • 


Corollary 

1.  ft(B)  is  always  non-singular  if  and  <d2  are  incommeasurable. 

2.  If  ui-|  and  u>2  are  both  rational  numbers,  we  can  make  ft( B)  nonsingular  by 
choosing 

B  <  (3.13) 

*'l 

where  T  is  the  period  denfined  in  (1.7). 

Proof.  Corollary  1  follows  directly  from  (3.12).  Corollary  2  is  proved  in 
Appendix  0. 

★ 

B.  Finding  the  initial  state  Xq 

Consider  the  1th  component  of  (3.5)  at  t  =  (2B+1)T^: 


xss.  ((2B+1  )Ti )  -  9if0(0)  +  ^  rgi,2k-l(0)  cos(2B+1>^Ti  +  9ij2k(0)  sin(2B+l)u)2T1j 

L  (3.14) 

Substituting  (3.11)  for  g^  ^(0)  in  (3.14),  we  obtain 

xssi  ((2B+1)Tt)  -  5T(B)gB_  =  5T(B)n-1(B)xss  (B) 

(3.15) 


r-1 


Xs.  (B)n'  (B)«(B) 

-ssi  - 


where 


(3.16) 


Since  [flT(B)]"‘  5(B)  in  (3.15)  remains  unchanged  for  all  1  =  1,2 . n,  we 

can  combine  all  "n"  components  of  x  (t)  from  (3.15)  Into  a  single  matrix  equation: 


xSSi«2B+l)T,) 

x,s2((2B*1>t,) 

xss  f(2B+l )T,) 


xss,“»  *ss,<T1>  Xas/2V 
xss2<«  xss2<Tl>  *ss2<2t!> 


xss  <°1  xss  <T1>  xss  <2T1> 
n  n  n 


^.((ffl+DT,) 
where  the  (2B+1)  -  vector 


X(B) 


Xss/2BT1> 

XSS2'2BT1> 


xssn(2BT,) 


80(B> 

B,(S) 

b2(b) 


b2B(b) 


B(B) 


(3.17) 


BCB)  a  [nT(B)]‘]  5(B) 


(3.18) 


depends  only  on  B. 

Observe  that  (3.17)  is  exact  provided  the  integer  bound  B  in  (3.4)  includes 

all  "M"  frequency  components  of  the  exact  steady-state  response  x$s(t)  in  (3.3). 

In  this  case,  the  entries  x$s(t),  t  =  0 Jj,2Tj ,. • . »2BT^  in  X(B)  can  be  obtained 

by  solving  (3.1)  using  xQ  as  the  initial  state. 

Since  xQ  is  precisely  what  we  are  seeking,  let  us  define  the  following 

system  of  "n"  related  equations  as  a  function  of  the  "n"  components 

xn  ,xn  ,...,xn  of  the  initial  vector 
U1  u2  un 


Jo  -  ^x0i  xq2  xq  =  x2^0^  xn(°)3T: 
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((2B+1  )Tj  ;xQ) 

x-|(0) 

X1 ^T1 *?0^ 

(2T  »  Xq )  •  •  • 

X1  *Xq) 

x2((2B+1)T1;x(J) 

x2(0) 

x2^T1 ;?0^ 

*2^]  »Xq)  •  •  • 

x2(2BTT;xQ) 

xn f(2B+l )T^ ;Xg) 

x„(0) 

xn^Tl '*0^ 

•  •  •  • 

•  •  • « 

xn^Tl  '~0^  '*• 

V281!^ 

x((2B+l )T] ;xQ) 


X(8;x0) 


e0(B) 
B-j  (B) 
82(B) 


B2B(B)J  JO 

* _ / 

8(B)  0 


F(xQ;B)  (3.19) 

where 

X^kTpXg),  k^O.l ,. . . ,2B+1  denote  the  1th  component  of  the  complete  solution 
x(t)  *  xtr(t)  +  xss (t )  of  (3.1)  starting  from  the  Initial  state  xQ,  for  t  *  0, 

Tj ,2Ti ,. . .2BT^ ,(2B+1 )T^ .  Observe  that  unlike  (3.17),  both,  x((2B+l ;x0)  and 
X(B;xn)  can  be  calculated  by  solving  (3.1)  numerically. 

*“  it  * 

Since  (3.19)  reduces  to  (3.18)  when  xQ  =  x^  ,  it  follows  that  Xq  can  be  found 
by  solving 


F(xQ;B)  -  0 

by  Newton- Raphson  iteration  as  in  (2.42),  with  F(-)  replaced  by  F(xn;B). 

The  Jacobian  matrix  Jp(XgJ/)  can  be  evaluated  by  the  sensitivity  network 
approach  [3]  if  the  circuit  is  given.  If  only  the  equation  (3.1)  is  given, 
the  Jacobian  matrix  can  be  calculated  from 


(3.20) 


(3.21) 


3x((2B+l )Tj ;Xg) 

wi*ere  ~ — - * —  can  be  evaluated  using  (2.55). 

*v0 

To  save  computation  time,  the  integer  bound  B  Is  chosen  to  include  only  the 
significant  frequency  components  In  (3.3).  In  this  case,  the  initial  state  xQ*(B) 
computed  from  (3.20)  will  depend  on  B  and  is  therefore  not  exactly  equal  to  xJf . 
Clearly, 
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(3.22) 


★  * 

* 

*o  tB>  -  2o 

as  B  -*•  B 

where  8*  denotes  the  Integer  which  is  large  enough  to  include  all.  frequency  com¬ 
ponents  of  (3.3). 

C.  Initialization  Guidelines 

Since  (3.1)  must  be  solved  many  times  numerically  from  t  *  0  to  t  =  (2B+1)T|. 
we  always  choose  oj-j  to  be  the  larger  of  the  2  input  frequencies.  Once  u-j  is 
identified,  we  can  calculate  F(xq-.B)  in  (3.19)  by  specifying  the  2  parameters 
B  and  Xq. 

a)  Choice  of  B. 

The  integer  bound  B  in  (3.5)  should  be  chosen  equal  to  at  least  the  order  of 
the  highest  significant  harmonics  of  w  In  the  steady-state  response.  It  is 
independent  of  u. .  This  Important  property  allows  us  to  analyze  a  large  class  of 
communication  circuits  where  the  "signal"  frequency  at  uu  is  much  smaller  than 
the  "carrier"  or  "pump"  frequency  u>i .  In  such  cases,  harmonics  of  ^  will  usually 
be  quite  small  even  though  the  input  signal  at  frequency  ui-|  is  usually  very  large 
(thereby  generating  many  higher  harmonics  of  co^)  so  that  accurate  answers  can  often 
be  obtained  with  B  3  3. 

If  the  order  of  the  highest  significant  harmonics  is  not  known,  we  simply 
make  an  intuitive  guess.  If  this  guess  is  unrealistic,  it  will  reveal  itself  in 
the  subsequent  error  analysis  (to  be  discussed  in  Section  III-D). 

b)  Choice  of  Xq. 

The  same  procedure  presented  in  Section  II-C  also  applies  here. 

D.  Termination  Guidelines. 

Recall  that  in  practice,  the  solution  xn  (B)  of  (3.20)  is  not  the  exact 
solution  Xq  .  Consequently,  we  must  val idate  this  answer  before  terminating. 

(1)  If  the  steady-state  response  x__(t)  Is  periodic  with  a  reasonably 

‘*>55  “  ""  ^ 

small  period  T  (see(1.7))  then  we  simply  solve  (3.1)  numerically  for  x(t;Xg  (B)) 

(with  Xq  (B)  as  initial  state)  and  verify  that  x(0;x0*(B))i  x(T;Xg*(B)). 

(2)  If  the  steady-state  response  x$s(t)  is  not  periodic,  or  if  it  is 
periodic  with  an  unreasonably  large  period  T,  we  can  estimate  the  error  with  the 
help  of  (2.51).  If  the  "approximate"  solution  x(t;Xg  (B))  Is  Indeed  close  to  the 
exact  solution  x(t;Xg  )  for  all  t  > 0,  then  it  follows  from  (2.43)  and  (2.51)  that 

lx(t;xQ*)  -  x(t;XQ*((B))l  <  l$(t)l  lxQ*  -  xQ*(B)l  (3.23) 

for  all  t  >  0. 

If  we  let  gk  (0)  and  9^(0)  denote  the  "exact"  (computed  using  (3.11)  with 
8  *  B  and  exact  Xq  )  and  "approximate"  (computed  using  (3.11)  with  approximate 
Xq  (B)}  values,  then 
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r 


[ 

j 

t 


B  B 

>5a*  -  ■  'Cso'dHI^on  -  Csoto)-j,52k-i((»3| 

i  '90*(01  -  90C0)I  *  kP  4_,(0)  -  g2k.,(0)l  *  jw  1^., (OH  (3.24) 

Since  the  coefficients  gk  (0)  and  9^(0)  can  be  interpreted  as  the  Fourier  co¬ 
efficients  associated  with- the  frequency  kwg  at  t  3  0,  it  is  reasonable  to  assume 

that  if  lgt(0)l  is  sufficiently  small  for  k  >  2B  +  1 ,  the  computed  initial  state 
*  * 

Xq  (B)  will  be  sufficiently  close  to  Xg  ,  and  hence 

g2k-l(0)  -§2k-l(0)’  k  *  (3.25) 


It  follows  from  (3.23),  (3.24),  and  (3.25)  that  we  can  approximate  (3.23)  by 


MtjxJ)  -  x(t;x*(B»l  <  l*(t)l{  ,g2k-1<0)1  +  c(B.x*(B))> 


(3.26) 


where  e(B,xQ*(B))  is  an  error  from  the  first  two  terms  in  (3.24). 

Even  though  the  right  hand  side  of  (3.26)  can  not  be  calculated  from  avail¬ 
able  data,  the  following  heuristic  procedure  has  been  used  successfully  in  all 
examples  we  have  investigated  so  far: 

(1)  Solve  (3.1)  numerically  for  x(t;x0*(8))  from  t  *  0  to  t  *  2(8+2)"^ . 

(2)  Calculate  gB  ,  1  *  l,2,...,n,  using  (3.9)  and  (3.11)  with  B  replaced 

by  8+2  and  with  x$s  (t)  replaced  by  xi(t;x0*(B))  for  t  3  0,T1,2T1 . 

2(B+2)T1 .  1 

(3)  If 

B+l£B+2  ^*52B+1^0^  +  1§2B+3^11 


is  smaller  than  some  prescribed  tolerance,  stop. 
Otherwise,  increase  B  and  start  all  over  again. 
E.  Summary  of  Two-Frequency  Algorithm 


Step  0.  Choose  to  be  the  larger  of  the  2  input  frequencies.  Specify  the 
integer  bound  B  (see  Section  III-C).  Set  j  ■  0. 

Step  1.  Choose  Initial  state  Xg  ■  x^  (for  j«0,  see  Section  II-Q)  and  solve 
(3.1)  numerically  to  obtain  x(T^),  x(2T^),  x(2BT^). 

Calculate  F(xgjB)  from  (3.19).  If  I F(xQ ;B)I  <  e,  where  e  Is  a 
sufficiently  small  positive  number,  call  xQ  »  Xg(B)  and  go  to  Step  4, 
Step  2.  Compute  the  Jacobian  matrix  Jp(Xg^)  in  (3.21)  (see  Section  II-C). 
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^Step  3.  Compute  Xg^+^  via  the  Newton-Raphson  iteration  (2.42)  with  F(*)  re¬ 
placed  by  F(Xg,B).  Go  to  Step  1  with  j  -*■  j+1. 

Step  4. 

Case  1 .  xs$(t)  is  periodic  with  small  period  T: 

Solve  (3.1)  for  x(t)  with  Xg  (B)  as  initial  state  from  t  *  0  to  t  *  T. 
Calculate 

ej  =  xj(°*?o*(8^  ‘  xj  j  m  1>2>*-'»n- 

Case  2.  xs$(t)  is  not  periodic  or  is  periodic  with  large  T: 

Solve  (3.1)  for  x(t)  with  Xg  (B)  as  initial  state  from  t  a  0  to 
t  =  2(8+2)^.  Calculate  g+ieg+2  as  defined  in  (3.27).  If 

max{e^,e2 . en^>e0  (easel) 

B+l eB+2  >  £0  (case  2^ 

where  Eg  is  a  sufficiently  small  positive  number.  Increase  B  and 
repeat  Steps  0-4. 

Step  5.  Stop. 


F.  Illustrative  Examples 

Example  1.  Puffing's  Equation  with  2  frequency  inputs: 

We  have  used  the  preceding  algorithm  to  solve  (2.64)  when  f(t)  contains  only 
2  Input  frequencies.  The  results  corresponding  to  3  different  combinations  of 
parameter  k,  c^,  c2>  and  f(t)  are  summarized  In  Table  5  for  B  *  9,  11,  13,  and 
17  respectively. 

Table  5.  Examples  Applying  the  Two-Frequency  Algorithm 


(1)  f(t)=0.5  cost 

(2)  f(t)*0.3  cost 

(3)  f(t)=(l+cos  0. 1 1 5t)  cost 

B 

+0.5 

cos  0. 81 t 

+1.5 

cos  0.1 15t 

1 

_ 

m 

V(B> 

BB 

tsa 

BB 

m 

9 

1.04898 

0.26642 

1.27285 

0.27251 

1.36899 

-0.34537 

1.11403 

0.64204 

1.24281 

0.17135 

1.39967 

0.00736 

III 

1.12986 

0.63906 

1.22548 

0.30906 

1.34835 

0.16875 

19 

1.11865 

0.63562 

1.21332 

0.33872 

1.35403 

0.15168 

16e17’ 

>0.0021 

16e17s 

>0.0083 

14e15=0'023 
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In  each  case,  an  error  estimate  using  (3.27)  is  calculated  and  the  results 
are  also  listed  in  Table  5.  For  example,  in  case  1,  we  have 

16e17  *  ,g32(0)l  +  ,?34C0)I 

£  lg32,l(°)i  +  l932,2^l  +  lg34,l^0^  +  lg34,2^0^ 

*  0.000853  +  0.000444  +  0.000244  +  0.000554  s  0.0021  (3.28) 


The  rate  of  convergence  for  these  3  cases  are  shown  in  Fig.  5.  The  conver¬ 
gence  rate  for  case  3  Is  not  as  good  as  the  cases  1  and  2  because  we  have  delib¬ 
erately  chosen  a  poorer  Initial  guess  for  constrast. 

The  steady-state  waveforms  corresponding  to  the  3  cases  listed  in  Table  5 
are  shown  in  Figs.  6(a),  7(a),  and  3(a),  respectively.  The  corresponding  frequency 
spectrum  calculated  by  the  FFT  algorithm  [17]  are  shown  in  Figs.  6(b),  7(b),  and 
8(b),  respectively.^  For  all  cases,  the  higher-order  harmonic  and  with  modulation 
components  are  negligible,  as  is  typical  in  many  practical  examples. 

Example  2.  Transistor  Modulator  Circuit: 

Consider  the  differential -pair  amplitude  modulator  circuit  shown  in  Fig.  9(a), 
where  e^t)  and  e2(t)  denotes  the  carrier  and  signal  input,  respectively.  Using 
the  algorithm  described  in  Appendix  A,  and  the  Ebers-Moll  circuit  model  [3]  shown 
In  Fig.  9(b)  for  the  transistors,  we  obtain  the  following  system  of  4  implicit 
differential-algebraic  equations  for  this  circuit: 


dvi  v 

•C  IT  ‘  Y  '  12  '  ls 


X(vrV„) 

[e  1  C  -1]  +  aIs[eXv4-l]  =  0 

\(v4+e,)  X(e,-V  )  .  A(v.-v  ) 

Ic[e  4  '  -1]  +  al  [e  1  cc  -1]  -  Ie[eXv4-l]  +  otl c[e  1  cc  -1] 


IsCe 


s—  "S 

X(v4+e2-VE) 


\v « 


1]  +  ctls[e  -1]  =  0 


(3.28) 


X(v4+e?-VF)  , 

-IsCeXv3-l]  +  als[e  4  i  E  -1]  +  yf-  [Vcc-VE+e2-v3]  *  0 

8 

-*■  ir  -  vi  * 0 


Note  that  the  first  3  equations  in  (3.28)  correspond  to  KCl  applied  at  nodes 
©,(£),  and®  respectively,  whereas  the  4th  equation  corresponds  to  KVL  applied 
around  the  loop  formed  by  the  Inductor  L. 

Note  that  unlike  In  the  Multi -Frequency  Algorithm,  the  generalized  Fourier  co¬ 
efficients  are  not  directly  ava. Table  in  this  algorithm. 
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Setting  e^t)  and  e2(t)  to  zero,  we  first  solve  (3.28)  for  the  following  dc 
operating  point: 

v.  «  0,  i,  =  0.155  x  10“3,  v,  =  0.259,  v.  -  0.241  (3.29) 

'Q  Q  Q 

Vie  then  choose  x  (0)  A  v,(0)  s  v,  and  x,(0)  A  i9(0)  =  i,  as  our  initial 

ini  1  —  1  'n  c  *“  c 

guess  x^  '  and  apply  the  two-frequencyyalgorithm  for  2  different  amplitudes  V2 

for  the  signal  egCt);  namely  V2  *  4.0  and  5.3,  respectively.  For  V2  =  4.0,  we 
choose  B  ■  6.  However,  for  V2  =  5.3,  we  choose  B  =  13  to  account  for  the 
additional  harmonics  that  are  likely  to  be  significant  in  view  of  the  larger  in¬ 
put  signal  amplitude.  In  both  cases,  our  algorithm  converges  in  2  iterations  and 
the  results  are  sumnarized  in  Table  6. 

Table  6.  Results  obtained  with  two-frequency  algorithm  using  a  4th  order  BDF 

-ft 

method  [3]  with  a  step  size  h  =  4u( 10  )  sec. 


Initial  State 

Error  Estimate 

case 

V'j  (0) 

i2(0) 

B+leB+2  for  V1 

B+l£ B+2  for  ^2 

C1) 

V2=4.0 

B=6 

-3.927 

0.2387(10~3) 

0.78O0"3) 

0 . 59(1 0~7 ) 

(2) 

V2=5.3 

B=13 

-3.422 

0.1138O0"3) 

0.23(10_1 ) 

0 . 95 ( 1 0~5 ) 

Using  the  2  initial  states  in  Table  6,  the  steady-state  waveforms  correspond¬ 
ing  to  the  modulator  output  voltage  Vg(t)  and  the  base-to-emitter-voltages  V^g(t)  for 
transistors  T-j  and  T3  are  shown  in  Figs.  10(a), (b),(c)  and  11  (a) ,(b) ,(c)*,  respec¬ 
tively.  Note  that  the  modulator  output  waveform  in  Fig.  10(a)  is  typical  of  those 
composed  of  a  carrier  and  2  side  band  frequencies  Even  the  waveforms 

VEfi(t)  in  Figs.  10(b)  and  (c)  are  quite  smooth,  indicating  the  absence  of  sig¬ 
nificant  higher-order  frequency  components.  Consequently,  very  accurate  results 
were  obtained  with  only  a  relatively  small  B  =  6. 

On  the  other  hand,  the  corresponding  waveforms  for  case  2  in  Fig.  11 
indicate  the  presence  of  many  more  frequency  components.  Consequently,  a  much 
larger  B  will  be  needed  to  obtain  results  of  acceptable  accuracy.  We  found 
B  ■  13  is  adequate  for  this  purpose. 

The  normalized  frequency  spectrum  corresponding  to  the  output  waveforms 
v q (t )  in  case  2  as  obtained  by  the  FFT  method  [17]  is  shown  in  Fig.  12. 
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IV.  Concluding  Remarks 

Two  efficient  algorithms  have  been  presented  for  finding  almost  periodic 
steady-state  response  of  nonlinear  circuits  and  systems. 

The  multi-frequency  algorithm  is  very  general  as  it  allows  any  number  of 
commensurable  or  incomnensurable  input  frequencies  . .  ,u>p.  Although  the 

output  normally  includes  only  harmonic  and  inter-modulation  frequency  components 
of  the  form  +  +  '"pi^p*  where  mjk  are  Integers,  other 

frequency  components,  such  as  subharmonics,  may  also  be  included  in  this  algorithm 
if  their  presence  are  suspected. 

The  two-frequency  algorithm  is  applicable  only  if  there  are  no  more  than  2 
input  frequencies.  This  restriction,  however,  is  more  than  compensated  by  its 
greatly  increased  computational  efficiency,  specially  when  the  steady-state 
response  contains  many  frequency  components.  That  this  algorithm  is  significantly 
better  than  algorithm  1  (when  applied  in  the  2-frequency  case)  is  best  seen  by 
comparing  the  number  of  respective  coefficients  in  Table  4.  Note  that  for 
B*10,  Algorithm  1  must  calculate  221  coefficients  whereas  Algorithm  2 
needs  to  calculate  only  21.  Note  that  2B  +  1  increases  only  by  2  as 
we  increase  B  by  1;  consequently,  the  two-frequency  algorithm  remains  comput¬ 
ationally  quite  efficient  even  with  a  larger  B,  thereby  allowing  stronger 
nonlinearities.  This  Is  particularly  useful  when  the  amplitude  of  the  higher- 
frequency  input  (c^  )  is  much  larger  than  that  of  u^,  as  is  coiimon  in  communication 
circuits  where  u-j  denotes  the  carrier  frequency  and  uig  denotes  the  signal 
frequency.  In  this  case,  the  number  of  significant  harmonic  components  of 
oi2  will  be  relatively  small  so  that  a  small  B  suffices. 

It  is  also  interesting  to  note  that  in  the  limiting  case  where  we  have  only 
one  frequency  input  (P-1),  then  (3.5)  reduces  to  x$s(t)  *  gQ(t) .  In  this  case, 
the  two-frequency  algorithm  reduces  to  the  usual  shooting  method  [3,11]. 

Certain  numerical  ill -conditioning  could  occur  in  the  Algorithm  1  when 
the  step  size  At  is  chosen  to  be  too  smal 1 .  The  ill-conditioning  problem  is  due 
to  loss  of  number  of  significant  digits  and  therfore  depends  strongly  on  the 
choice  of  the  computer. 

Finally  we  remark  that  if  the  steady  state-solution  is  not  periodic  so  that 
the  brute-force  method  is  impractical  (since  we  must  in  theory  integrate  for  alj_ 
t  >.  0),  or  if  the  nonlinearity  is  not  sufficiently  weak  for  the  Perturbation  and 
Vol terra  series  methods  to  converge,  then  our  algorithms  are  presently  the  only 
methods  available  for  finding  steady-state  solutions,  let  alone  their  good  compu¬ 
tational  efficiency. 
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V.  Appendix 


Appendix  A.  Explicit  Formula  for  Reduced  System  of  Implicit  Equations 

Let  N  be  a  nonlinear  network  containing  voltage  or  current-controlled 
2-terminal  resistors,  voltage-controlled  2-terminal  capacitors,  current- 
controlled  2-terminal  inductors,  as  well  as  independent  and  controlled  sources. 
Mutual  couplings  are  allowed  so  long  as  they  are  restricted  to  elements  belonging 
to  the  same  class.  Let  each  independent  source  be  considered  as  part  of  a 
"composite"  branch  as  in  [3].  Adoping  the  notations  in  Section  17-2  of  [3],  we 
obtain  the  following  tableau  equation  for  N: 


n-1  KCL  equations  { 

A 

0  C 

i 

b  KVL  equations  { 

0 

1  -at 

V 

~ 

- 

b  elements  consitu-{ 
tive  relations 

~v 

A 

AJ 

E 


(A-l ) 


Equation  (A-l)  consists  of  a  system  of  (n-1)  +  2b  implicit  equations  of  the  form 
(1.1)  where  "b"  denotes  the  number  of  composite  branches  and  "n"  denotes  the 
number  of  nodes.  Our  goal  in  this  section  is  to  derive  an  equivalent  system  of 
implicit  equations  containing  fewer  number  of  equations  and  variables  for  an 
important  subclass  of  networks. 

In  particular,  we  assume  that  N  contains  no  loops  of  capacitors  and  indepen¬ 
dent  voltage  sources,  no  cut  sets  of  inductors  and  independent  current  sources 
and  that  all  controlled  sources  are  current  sources  depending  on  either  resistor 
or  capacitor  voltages.  Consequently,  there  always  exists  a  normal  tree  T  con¬ 
taining  all  capacitors  and  no  inductors  [3]. 

If  we  let  ig  and  tz  denote  the  current  and  voltage  vectors  of  al_l_  inductors 
in  N,  and  let  i^  and  y-j  denote  the  current  and  voltage  vectors  of  the  remaining 
elements,  then  (A-l)  can  be  recast  as  follows: 


(A-2) 


where  the  reduced  incidence  matrix  A  is  similarly  partitioned  into  A  =  [A-j  A£], 
and  where  Mi^)  denotes  the  Incremental  inductance  matrix.  Substituting 
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into  the  last  equation  in  (A-2),  we  obtain  the  following  reduced  system  of 
equations: 


tzlz  +  ^llb^n  =  illb!l  •  ^l?(~V~v)  +  " 
*2  =  &  +  e2  =  l(i2)i2 


(A-4) 

(A-5) 


Let  vT  denote  the  branch  voltage  vector  associated  with  the  normal  tree  T, 
and  let  v^  denote  the  corresponding  cotree  voltages.  Since  all  capacitors  are 
assigned  in  T,  v^  is  a  subvector  of  Vj.  Similarly,  since  all  inductors  are 
assigned  in  the  cotree,  v2  is  a  subvector  of  v^.  Let  the  reduced  incidence 
matrix  A  be  partitioned  accordingly  into  Aj  and  A^,  so  that  KVL  assumes  the  form 

(A-6) 


Since  the  columns  of  kj  correspond  to  tree  branches,  Aj  is  non-singular  [3]. 
Hence  we  can  solve  for  the  node- to-da turn  voltage  vector  vn  from  (A-6)  to  obtain 

\  Evr  -  Er]  (A-7) 

vL  -  ?Mr’rvT  -  y  *  el  (a-s) 

Substituting  (A-7)  and  (A-8)  into  (A-4)  and  (A-5),  and  denoting  the  inductor 
current  vector  i2  by  i^,  we  obtain 

(A-9) 
(A-10) 

where 

g(vc,vr)  a  g(y  .v)| 

|y  =  [vT  vLy 

and  Vj  is  given  by  (A-8), 

Equations  (A-9)-(A-lQ)  constitute  a  reduced  system  of  Implicit  equations  in 
terms  of  the  state  variables  x  A  [v^i^]T  and  the  non-state  variables  contained 
within  Vy. 


A21L  +  ^l~b~l^~T^  ^  -p)  +  ~1  Xb~l  +  -li^c’-T^  ~  ~~  ~  ~ 

A2[Aj  ]  _1 (Vj  -  et)  +  E2  -  L(iL)iL  =  g 


Equation  (A-9)  can  be  interpreted  as  the  nodal  equation  of  N  with  all 
inductor  currents  i^  considered  as  independent  sources,  and  with  all  node- 
to-datum  voltages  expressed  in  terms  of  the  normal  tree  voltage  vector  Vj. 
Similarly,  (A-10)  can  be  interpreted  as  the  fundamental  loop  equations  (relative 
to  the  normal  tree  T)  formed  by  the  inductor  links.  These  interpretations  allow 
us  to  write  down  the  reduced  system  of  implicit  equations  of  simple  nonlinear 
networks  --  such  as  that  considered  in  section  III  —  by  inspection.  If  N 
contains  loops  of  capacitors  and  independent  voltage  sources,  or  cut  sets  of 
inductors  and  independent  current  sources,  the  above  procedure  can  be  easily 
generalized  by  first  expressing  the  cotree  capacitor  voltages  in  terms  of  tree 
capacitor  voltages,  and  tree  inductor  currents  in  terms  of  cotree  inductor 
currents  [3]. 

APPENDIX  B 

The  matrix  rTr  in  Theorem  3  is  non  singular  if,  and  only  if, 
a*  4  2nrr 

TvvT 

where  V  is  a  Zx(2N+l)  matrix,  where  Z  >  2N+1. 

Proof.  rTf  is  clearly  symnetric  and  positive  semi-definite  and  hence  its 
eigenvalues  Xx,  X2,...,Xn  are  real  and  non-negative.  Hence  det  (r  r) 

»  Xj,  X2,...Xq  +  0  all  eigenvalues  are  positive  •  rTr  is  positive  definite  • 
xT(CTr)  x  >0  for  all  x  t  0 
Hence,  we  have 

det  (£T£)  t  0  **  columns  of  E  are  linearly  independent.  (B-2) 

The  matrix  E  in  (2.17)*  can  be  recast  as  follows: 


Vs  1  g-j29-j  ^ 


1  1  1  0  0  ...  0  0 

ej®N  e"j®N  o  [o.5  -j0.§]  ...  0  0 

.  ej20N  e“J’20N  0  0.5  j0.5  ...  0  0 


1  ejZ9l  e-jZ9l  .  .  .  ejZ0N  e’jZ0N 


0.5  -j0.5 
0.5  j0.5 


where  9^  ^  ->.At.  Since  0  is  clearly  non-singular, 
Columns  of  r  are  linearly  independent 


Columns  of  £'  are  linearly  independent 


e.  ±  0k  i2mr 


e"J  i  t  1  for  any  i  f  k 
9^  t  ±2mr,  for  any  integer  n 


•  0.  ±  8^  f  ±2nir  for  any  integer  n 

It  follows  from  (8.2)  and  (B.4)  that 

det  (rTr)  t  0  •  Vi  it  i  \).  at  ^  ±2mr  «•>  bX  f  -r—S 
'  -  I  N  |v  j  sv 


APPENDIX  C 

The  matrix  0(8)  in  (3.9)  is  non  singular  if,  and  only  if,  there  exists  an  in¬ 
teger  l_2  such  that 

7T  *  r-  •  L,  -  1,2,. ..,28  (C-l ) 

U1  L1  1 

Proof .  We  can  recast  0(B)  in  (3.9)  as  follows: 


1  1 


0(B)  = 


1  eJ2<p  fi-j2iJ> 


ej84>  e-jB* 
e j 2B4>  e-j2Bd> 


0.5  -j0.5 
0.5  j0.5 


1  ei2B<P  e'j2B4> 


5j2B20  e-j2B2<p 


0.5  -j0.5 
0.5  j0.5 


(pTB) 

where  $£<*>2Tv  ^ince  8  Is  non-singular, 

0(B)  is  non-singular 

~  Columns  of  O' (8)  are  linearly  independent 
«*  e^k<1>  -  e^  t  0  for  1,k  a  0,±1,  ±2,...tB,  i  f  k 

„  1  .  ejLl<J>  f  Q  for  l,  «  1,2,...  ,2B 

«*  L1  $  f  2L2jt  for  any  integer  L2 

Substituting  $  a  u^T-j  *  2ttu)2/u>-j  into  (C-3),  we  obtain: 

0(B)  is  non-singular 


(C-2) 


(C-3) 


*  Lj  (u^/oi^ )  t  L,. 


A-4 


APPENDIX  D 

If  and  l>2  are  rational  numbers,  we  can  make  ft(B)  non- singular  by  choosing 


Proof.  Let  a  m^/n^  and  w2  A  m2/n2  be  irreducible  fractions.  Then 
*  2ir{n.|/n^ ),  T2  *  27r(n2/m2),  and  T  =  27r(n/m),  where  n  «  L.C.M.  {n1  ,n2>  and 
m  *  G.C.D.{m.j,m2}.  Hence  m-,  »  mm.,'  ,  m2  =  mm2’  ,  n1  *  kny  ,  n2  «  kn2'  ,  and 
n  =  kn^'  n2'  for  some  integer  k. 

Now,  since  G.C.O.^'  ,m2’  }  *  1,  G.C.D.&n^  .n^  }  =  1,  and  G.C.D.{m2'  ,n2'  }*  1 
we  have  G.C.D.Im^  n^ , rn^.n^  *  1.  It  follows  that  if 

max{L|}  =  2B  <  m^  n2'  (D«2) 


’Q)  *  ■  l'(S)  *  L,(?v) ' ,nte3er 

But  (D.2)  is  equivalent  to 

/m,‘  m\  /n,'  n.\ 

m1  n2  Tl  *  <"1'n2l(n1/m1)(m/n)T  =  \-l—jLL-±\  j  , 


CD-3) 


2BT-J  <  mi  'n2'Ti  *  m-j 1  n2 1  (ni/m-j )  (m/n)T 


Hence,  if  28T.|  <  T,  then  (C.l)  holds. 


n,'  kn,1 

T  *  k^TaT  (D*4) 
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FIGURE  CAPTIONS 


-J 

Fig.  1  Steady-state  waveform  for  Duffing’s  equation  x  +  O.lx  +  2x  +  x 

-  0.4  cos  t  +  0.4  cos  Q.35t  +  0.4  cos  0.1 55t. 

(a)  Solid  waveform  represents  approximate  solution  obtained  with 

multi -frequency  algorithm.  Solid  dots  denote  solution  obtained  by 

numerical  solution  of  (2.64)  starting  from  x  *  =  (0.69667,-0.18304) 

*•0 

from  Case  1. 

(b)  Discrete  frequency  spectrum  obtained  from  (2.27). 

Fig.  2  Steady-state  waveform  for  Duffing' s  equation  x  +  O.lx  +  2x  +  x 
»  0.4  cos  t  +  +  0.4  cos  0.85t  +  0.4  cos  0.17t. 

(a)  Solid  waveform  represents  approximate  solution  obtained  with 
multi -frequency  algorithm.  Solid  dots  denote  solution  obtained  by 
numerical  solution  of  (2.64)  starting  from  xQ*  =  (0.79298,-0.13834) 
from  Case  2. 

(b)  Discrete  frequency  spectrum  obtained  from  (2.27). 

Fig.  3  The  rate  of  convergence  for  Case  1  (shown  dotted)  and  Case  2  (shown 
solid)  in  the  Examples  in  Table  2.  Horizontal  axis  indicates  the 
Iteration  number  j.  Vertical  axis  indicates  the  error  computed 
at  the  jth  Iteration  using  (2.67). 

Fig.  4  Geometrical  interpretation  of  (3.4)  for  B  *  1,2,. ...10.  Each  solid 

dot  denotes  one  frequency  component  m1  ^  +  m2ku>2. 

Fig.  5  The  rate  of  convergence  for  Cases  1,  2,  and  3  In  the  Examples  in 

Table  3.  Horizontal  axis  Indicates  the  Iteration  number  j.  Ver¬ 
tical  axis  indicates  the  error  estimated  by 

e(j)  „ 

Fig.  6  (a)  Steady-state  waveform  for  Duffing' s  equation  x  +  0.06x  +  x  +  x3 

■  0.5  cos  t  +  0.5  cos  0.81  t  (Case  1  ) 

(b)  Normalized  frequency  spectrum  of  (a). 

Fig.  7  (a)  Steady-state  waveform  for  Duffing's  equation  x  +  0.05x  ♦  x  +  x3 

■  0.3  cos  t  +  1.5  cos  0.115  t  (Case  2) 


(b)  Normalized  frequency  spectrum  of  (a). 

Fig.  8  (a)  Steady-state  waveform  for  Duffing' s  equation  x  +  O.lx  +  x  +  x 

*  (1  +  cos  0.115  t)  cos  t  (Case  3) 

(b)  Normalized  frequency  spectrum  of  (a). 

Fig.  9  (a)  Differential-pair  amplitude  modulator  circuit.  Vcc  *  10V, 

VE  *  5V,  L  -  2  mH,  C  «  500  pF,  R,_  =  20  kQ,  Rg  *  15  kfl,  e^t)  » 

*  0.01  cos  0. 115(1 0®)t  and  e2(t)  =  V2  cos  0.115(106)t. 

(b)  Ebers-Moll  ..transistor  circuit  model  with  the  2  diodes  described 
by  *dk  "  !sCe  “  U,  Is  «  10  SA,  X  *  40,  a  *  0.99. 

Fig.  10  (a)  Steady-state  output  voltage  waveform  vQ(t)  for  Case  1 :  carrier 

signal  e^(t)  *  0.1  cos  10®  t,  input  signal  e2(t)  *  4.0  cos  0.115(l0®)t 

(b)  Corresponding  base-to-emitter  voltage  waveform  vEg(t)  for 
transistor  T 

(c)  Corresponding  base-to-emitter  voltage  waveform  vEg(t)  for 
transistor  Tg. 

Fig.  11  (a)  Steady-state  output  voltage  waveform  VQ(t)  for  Case  2:  carrier 

signal  e^t)  -  0.1  cos  10®  t,  input  signal  eg(t)  =  5.3  cos  0.115(10®)t. 

(b)  Corresponding  base-to-emitter  voltage  waveform  vgg(t)  for 
transistor  T-j . 

(c)  Corresponding  base-to-emitter  voltage  waveform  vEg(t)  for 
transistor  Tg. 

Fig.  12  Normalized  frequency  spectrum  for  the  modulator  output  voltage  wave¬ 
form  in  Fig.  11(a)  (Case  2) 


ERROR  *  € 


l 


1 


-10  -9  -8  -7  -6  -5  -4  -3  -2  -I  0  I  2  3  4  5  6  7  8  9  10 


Fig.  4 


Normalized  Frequency  Spectrum  of  x(t) 


Normalized  Frequency  Spectrum  of  x(t) 


Normalized  Frequency  Spectrum  of  x(t) 
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